# Copyright (C) 2010 Zajcev Evgeny <zevlg@yandex.ru>

# C-fast implementation of ``is_inside`` using winding number method +
# grid optimization technique

# Take a look at testing suite at ``tests/test_is_inside.py`` as well

import operator
from itertools import takewhile, chain
from functools import partial

DEF DBLMAX = 1000000.0

ctypedef struct vertex:
    double x
    double y

ctypedef struct contour:
    int n
    vertex* vertices

ctypedef struct polygon:
    int n
    contour* contours
    int* holes

ctypedef struct boundbox:
    double xmin
    double ymin
    double xmax
    double ymax

cdef enum: GB_UNKNOWN, GB_INSIDE, GB_OUTSIDE
ctypedef struct gridbox:
    int type
    boundbox bounds

# Winding number algorithm
cdef inline int is_inside_contour(contour* cnt, double x, double y):
    cdef int i, j = cnt.n-1, wn = 0
    cdef vertex *vi, *vj
    for i in range(cnt.n):
        vi = cnt.vertices+i; vj = cnt.vertices+j
        if vj.y <= y:
            if vi.y > y and (vi.x-vj.x)*(y-vj.y) > (x-vj.x)*(vi.y-vj.y):
                wn += 1
        elif vi.y <= y and (vi.x-vj.x)*(y-vj.y) < (x-vj.x)*(vi.y-vj.y):
            wn -= 1
        j = i
    return wn

def bymin(a,b): return map(min, a, b)
def bymax(a,b): return map(max, a, b)

DEF GRIDSIZE = 40

cdef class Polygon:
    cdef:
        polygon poly
        boundbox pbounds
        gridbox* grid                   # optimization grid
        double gxstep, gystep           # steps in grid
        public object contours          # P-contours
        public ghits, wnhits

    def __cinit__(self):
        self.poly.n = 0
        self.poly.holes = self.poly.contours = NULL
        self.contours = []

        self.grid = <gridbox*>malloc(GRIDSIZE*GRIDSIZE*sizeof(gridbox))

        cdef boundbox* sbb = &self.pbounds
        sbb.xmin = sbb.ymin = DBLMAX
        sbb.xmax = sbb.ymax = -DBLMAX

        self.ghits = self.wnhits = 0

    def __dealloc__(self):
        cdef polygon* sp = &self.poly
        free(self.grid)
        for i in range(sp.n):
            free(sp.contours[i].vertices)
        if sp.contours: free(sp.contours)
        if sp.holes: free(sp.holes)

    def add_contour(self, cnt):
        """Append a contour CNT to the Polygon.
        If HOLE is True than given contour denotes a hole.
        CNT is a list of points, point is 2-place sequence."""
        self.contours.append(cnt)

    def optimize(self):
        """Optimize the Polygon.
        Highly recommended if heavily using ``is_inside``."""
        cdef polygon* sp = &self.poly
        sp.n = len(self.contours)
        sp.contours = <contour*>malloc(sizeof(contour)*sp.n)
        sp.holes = <int*>malloc(sizeof(int)*sp.n)

        # Create contours and bound box
        cdef boundbox* bb = &self.pbounds
        cdef vertex* vp
        cdef int icnt = 0
        for cnt, hole in self.contours:
            vn = len(cnt)
            vp = <vertex*>malloc(sizeof(vertex)*vn)
            for i in range(vn):
                vp[i].x, vp[i].y = cnt[i]
            sp.contours[icnt].n = vn
            sp.contours[icnt].vertices = vp
            sp.holes[icnt] = hole
            icnt += 1

            # Adjust global bound box
            # NOTE: This bound box is rough because holes may affect
            # boundaries and we does not handle them
            if hole: continue
            bb.xmin, bb.ymin = reduce(bymin, cnt, (bb.xmin, bb.ymin))
            bb.xmax, bb.ymax = reduce(bymax, cnt, (bb.xmax, bb.ymax))

        # Create an optimization grid
        cdef gridbox* gbp
        cdef boundbox* gbb
        cdef vertex gbverts[4]
        cdef contour gcnt
        gcnt.n = 4; gcnt.vertices = gbverts

        self.gxstep = (bb.xmax - bb.xmin)/GRIDSIZE
        self.gystep = (bb.ymax - bb.ymin)/GRIDSIZE
        for yo in range(GRIDSIZE):
            for xo in range(GRIDSIZE):
                gbp = &self.grid[yo*GRIDSIZE+xo]
                gbb = &gbp.bounds
                gbb.xmin = bb.xmin+xo*self.gxstep
                gbb.xmax = gbb.xmin + self.gxstep
                gbb.ymin = bb.ymin+yo*self.gystep
                gbb.ymax = gbb.ymin + self.gystep

                vxs = [(gbb.xmin, gbb.ymin), (gbb.xmin, gbb.ymax), \
                       (gbb.xmax, gbb.ymin), (gbb.xmax, gbb.ymax)]
                vxsi = map(self._pis_inside, vxs)

                # Check that none of polygon vertex lies inside
                # gridbox, otherwise gridbox type is UNKNOWN
                for i in range(len(vxs)):
                    gbverts[i].x, gbverts[i].y = vxs[i]
                pp_in_gcnt = False
                for cnt, _ in self.contours:
                    for cp in cnt:
                        if is_inside_contour(&gcnt, cp[0], cp[1]):
                            pp_in_gcnt = True
                            break
                    if pp_in_gcnt: break

                if pp_in_gcnt: gbp.type = GB_UNKNOWN
                elif all(vxsi):
                    gbp.type = GB_INSIDE
                elif all(map(operator.not_, vxsi)):
                    gbp.type = GB_OUTSIDE
                else:
                    gbp.type = GB_UNKNOWN

        return self

    def bounds(self):
        """Return overall bounds of the Polygon."""
        cdef boundbox* sbb = &self.pbounds
        return sbb.ymin, sbb.xmin, sbb.ymax, sbb.xmax

    def _pis_inside(self, pnt):
        return self._is_inside(pnt[0], pnt[1])

    cdef inline int _is_inside(self, double x, double y):
        cdef int r = False
        for i in range(self.poly.n):
            if is_inside_contour(&self.poly.contours[i], x, y):
                r = not self.poly.holes[i]
        return r

    cpdef bint is_inside(self, double x, double y):
        """Return True if point X,Y is inside Polygon.
        Contours checked in order of appendence."""
        cdef boundbox* bb = &self.pbounds

        # Check global bounds first
        if x < bb.xmin or x > bb.xmax \
               or y < bb.ymin or y > bb.ymax:
            return False

        # Check coressponding gridbox
        cdef int yoff = <int>((y - bb.ymin)/self.gystep)
        cdef int xoff = <int>((x - bb.xmin)/self.gxstep)
        cdef gridbox* gb = &self.grid[yoff*GRIDSIZE+xoff]
        if gb.type is not GB_UNKNOWN:
            self.ghits += 1
            return gb.type == GB_INSIDE

        # Check the contours
        self.wnhits += 1
        return self._is_inside(x, y)
